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ABSTRACT 

Angular momentum transport owing to hydrodynamic turbulent convection is studied using local three di- 
mensional numerical simulations employing the shearing box approximation. We determine the turbulent vis- 
cosity from non-rotating runs over a range of values of the shear parameter and use a simple analytical model 
in order to extract the non-diffusive contribution (A-effect) to the stress in runs where rotation is included. Our 
results suggest that the turbulent viscosity is of the order of the mixing length estimate and weakly affected by 
rotation. The A-effect is non-zero and a factor of 2-4 smaller than the turbulent viscosity in the slow rotation 
regime. We demonstrate that for Keplerian shear, the angular momentum transport can change sign and be 
outward when the rotation period is greater than the turnover time, i.e. when the Coriolis number is below 
unity. This result seems to be relatively independent of the value of the Rayleigh number. 
Subject headings: accretion, accretion disks - convection - stars: rotation - Sun: rotation - turbulence 



1. INTRODUCTION 

Turbulence due to the convective instability is thought to 
account for much of the angular momentum transport in the 
outer layers of the Sun and other stars with convection zones 
(e.g. Riidiger 1989; Rudiger & Hollerbach 2004). In the pres- 
ence of turbulence the fluid mixes efficiently and diffusion 
processes occur much faster than in its absence. This effect 
is usually parameterized by a turbulent viscosity v t that is 
much larger than the molecular viscosity v. Often the value 
of v t is estimated using simple mixing length arguments with 
Vt = u rms l/3, where M lms is the rms velocity of the turbulence 
and I = omltH where qtmlt is a parameter of the order unity 
and H is the vertical pressure scale height. Numerical results 
from simpler fully periodic isotropically forced systems sug- 
gest that the mixing length estimate gives the correct order 
of magnitude of turbulent viscosity (e.g. Yousef et al. 2003; 
Kapyla et al. 2009a; Snellman et al. 2009). However, it is 
important to compute v t from convection simulations in order 
to see whether the results of the simpler systems carry over 
to convection. Furthermore, it is of interest to study whether 
the small-scale turbulent transport can be understood in the 
light of simple analytical closure models that can be used in 
subgrid-scale modeling. Measuring v t and its relation to aver- 
aged quantities, such as correlations of turbulent velocities, is 
one of the main purposes of our study. 

In addition to enhanced viscosity, turbulence can also lead 
to non-diffusive transport. The a-effect (e.g. Rrause & Radler 
1980), responsible for the generation of large-scale magnetic 
fields by helical turbulence, is one of the most well-known 
non-diffusive effects of turbulence. An analogous effect ex- 
ists in the hydrodynamical regime and is known as the A- 
effect (Rrause & Rudiger 1974). The A-effect is proportional 
to the local angular velocity and occurs if the turbulence is 
anisotropic in the plane perpendicular to the rotation vector 
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(Riidiger 1989). The existence of the A-effect has been estab- 
lished numerically from convection simulations (e.g. Pulkki- 
nen et al. 1993; Chan 2001; Kapyla et al. 2004; Rudiger et 
al. 2005) and simpler homogeneous systems (Kapyla & Bran- 
denburg 2008). 

If, however, both shear and rotation are present it is difficult 
to disentangle the diffusive and non-diffusive contributions. 
This is particularly important in the case of accretion disks 
where the sign of the stress determines whether angular mo- 
mentum is transported inward or outward. Convection is com- 
monly not considered as a viable angular momentum transport 
mechanism in accretion disks since several studies have indi- 
cated that the transport owing to convection occurs inward 
(e.g. Cabot & Pollack 1992; Ryu & Goodman 1992; Stone & 
Balbus 1996; Cabot 1996; Rudiger et al. 2002). Furthermore, 
in an influential paper, Stone & Balbus (1996, hereafter SB96) 
presented numerical simulations of hydrodynamic convection 
where the transport was indeed found to be small and directed 
inward on average. This result was used to provide additional 
evidence for the importance of the magneto-rotational insta- 
bility (Balbus & Hawley 1991) as the main mechanism pro- 
viding angular momentum transport in accretion disks. Al- 
though we agree with the conclusion that hydrodynamic tur- 
bulence is ineffective in providing angular momentum trans- 
port, we are concerned about the generality of the result of 
SB96. There are now some indications that hydrodynamic 
turbulence may not always transport angular momentum in- 
ward (cf. Lesur & Ogilvie 2010). 

In order to approach the problem from a more general per- 
spective, Snellman et al. (2009) studied isotropically forced 
turbulence under the influence of shear and rotation and found 
that the total stress, corresponding to the radial angular mo- 
mentum transport in an accretion disk, can change sign as ro- 
tation and shear of the system are varied in such a way that 
their ratio remains constant. They found that the stress is pos- 
itive for small Coriolis numbers, corresponding to slow rota- 
tion. In what follows we show that outward transport is possi- 
ble also for convection in a certain range of Coriolis numbers. 
When rotation is slow, the Reynolds stress is positive, corre- 
sponding to outward transport. In the regime of large Coriolis 
numbers the Reynolds stress changes sign and is directed in- 
ward, as in the study of SB96. 

The remainder of the paper is organized as follows: our 
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numerical model is described in Section [2] and the relevant 
mean-field description in Section [3] The results and the con- 
clusions are given in Sections[4]and[5] respectively. 

2. MODEL 

Our model is the same as that used in Kapyla et al. (2008), 
but without magnetic fields. We solve the following set of 
equations for compressible hydrodynamics, 
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where D/Dt = d/dt + ( U+ Uq) ■ V is the total advective deriva- 
tive, Uq = (0, S x, 0) is the imposed large-scale shear flow, v 
is the kinematic viscosity, K is the heat conductivity, p is the 
density, U is the velocity, g is the gravitational acceleration, 
and = Qoz is the rotation vector. The fluid obeys an ideal 
gas law p = pe(y - 1), where p and e are pressure and internal 
energy, respectively, and y - c P /cy — 5/3 is the ratio of spe- 
cific heats at constant pressure and volume, respectively. The 
specific internal energy per unit mass is related to the tem- 
perature via e = cyT. The traceless rate-of-strain tensor S is 
given by 

s. ! -i«//, ; + /:, / ,-U ; v-r. (4) 

where commas denote partial differentiation. 

We use a three layer, piecewise polytropic stratification 
with constant gravity g = -gz. The vertical variation of g 
in real accretion disks is neglected in our local model. The 
positions of the bottom of the box, bottom and top of the 
convectively unstable layer, and the top of the box, respec- 
tively, are given by (Z1.Z2.Z3.Z4) = (-0.85,0, 1, \ .\5)d, where 
d is the depth of the convectively unstable layer. Initially the 
stratification is piecewise polytropic with polytropic indices 
(m\,m2,mi,) = (3, 1, 1). The horizontal extent of the box is 
twice the vertical extent, i.e. L x /d = L y /d = 2L z /d = 4. 
The thermal conductivity has a vertical profile that maintains 
a convectively unstable layer above a stable layer at the bot- 
tom of the domain. A prescribed temperature gradient at the 
base maintains a constant heat flux into the domain. The last 
term in Equation 0} describes an externally applied cooling 
with 

^ T coo l(z)' ^ ^ 

where eo is the internal energy at z\ and t coo i(z) is a cooling 
time which is constant for z > Z3 and smoothly connects to 
the lower layer where r coo i — > 00. 

We use impenetrable stress-free boundary conditions at the 
top and bottom boundaries for the velocity 
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The temperature gradient at the bottom of the domain is given 
by 

dT = -g/Cy 

dz (mi + l)(y-l) 

where ni\ — 3 is the polytropic index at z\. All quantities are 
periodic in the y-direction, whereas shearing periodic condi- 
tions (Wisdom & Tremaine 1988) are used in the x-direction. 



The same setup has been used in earlier work to model con- 
vection in local patches in a star, but here we apply it also 
to a layer near the surface of an accretion disk. The source 
of heating in the midplane of the disk is not specified and 
is instead assumed given. This is appropriate for addressing 
the more general question about the direction of angular mo- 
mentum transport once there is convection in the absence of a 
magnetic field. 

We employ the Pence. Code 5 , which is a high-order finite 
difference code for solving the compressible MHD-equations. 
The bulk of our simulations were performed at a moderate 
resolution of 128 3 grid points. In a few cases we study the 
behavior of the solutions at higher resolutions (up to 1024 3 ), 
see Figure [TJ for a snapshot of a high resolution simulation, 
showing the vertical velocity on the periphery of the domain. 



2.1. Dimensionless quantities and parameters 
Non-dimensional quantities are derived such that 
d = g=p = c v = l, 



(8) 



where po is the initial density at zi- The units of length, time, 
velocity, density, and entropy are 

[x] = d, [t] = yfdfi, [17] = y[d~g, [p] = po, [s] = c P . (9) 

The system is characterized by several non-dimensional pa- 
rameters. We define the Prandtl number and the Rayleigh 
number as 

ft.i, Ra=^(-i£), (10, 

where xo = K/ (p m cp) is a reference value for the thermal dif- 
fusivity, and p m is the density at the mid layer z m = |fo - zi). 
The entropy gradient, measured at z m in the non-convecting 
hydrostatic state, is given by 
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where V - V a( j is the superadiabatic temperature gradient with 
V a d = 1 — 1/y, V = (dln77<91np)j. m , and Hp is the pressure 
scale height at z m . The amount of stratification is determined 
by the parameter - (j - l)eo/(gd), which is the pressure 
scale height at the top of the domain normalized by the depth 
of the unstable layer. We use in all cases £0 = 5. which results 
in a total density contrast of about 23. The Mach number in 
our simulations is of the order of 0.1 or less. 

We define the Coriolis and shear numbers, describing the 
strengths of rotation and shear, respectively, as 
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2Q 



Sh = 



(12) 



where u rms is the rms value of the turbulent velocity averaged 
over the full volume, kf = 2n/d is an estimate of the energy 
carrying wavenumber, and d is the depth of the convectively 
unstable layer. The Reynolds number is given by 

^rms 

vkf 
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For accretion disk applications it is convenient to define also 
the relative shear rate, 



S_ 



(14) 
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Fig. I. — Snapshot of the vertical velocity from run D5 with Re as 648. The sides of the box show the periphery of the domain whereas the top and bottom slices 
show U z at vertical heights z. = 0.95d and z. = 0, respectively. Here, Co = 0.26 and Sh = -0.20. See also http : //www . helsinki . £i/~kapyla/movies . html 
for an animation. 



which describes the rotation profile of the disk, where the lo- 
cal angular velocity varies like £2 oc r~ q . For Keplerian disks 
we have q — 3/2. 

The majority of the simulations ran for 10 3 time units 
(Af = 10 3 y/d/g), which corresponds to roughly 250 convec- 
tive turnover times; r tuln = (Wrms&f) • For the analysis the first 
50 turnover times of the simulations were usually discarded in 
order to minimize the effects of initial transients. The highest 
resolution runs with 512 3 and 1024 3 grid points were started 
by remeshing from snapshots of lower resolution runs and ran 
for 60 and 14 turnover times, respectively. 

Error bars are estimated by dividing the time series into 
three equally long parts and computing averages for each part 
individually. The largest departure from the average over the 
full time series is taken to represent the error. 

3. MEAN-FIELD INTERPRETATION 

In mean-field hydrodynamics the velocity field is decom- 
posed into its mean and fluctuating parts, 

U=U + u, (15) 

where the overbar denotes averaging and lowercase u the fluc- 
tuation. In the present paper we consider horizontal averaging 
so that the mean quantities depend only on z. We define the 
Reynolds stress as 

Rij = UiUj. (16) 

Fluctuations of the density are here ignored for simplicity. 

The Reynolds stress is often described in terms of the 
Boussinesq ansatz which relates the stress to the symmetrized 



TABLE 1 

Summary of the Runs with Variable q. Here, Ma = Krms/G?^) AND 



R\y = {Rjcy)/W%ms! WHERE (R„) IS THE VOLUME AVERAGE OF THE STRESS OVER THE 
CONVECTIVELY UNSTABLE LAYER. Pr = 1 AND Ra = 10 S IN ALL RUNS. 



Run q Co Sh Re Ma grid 



AO - 0.00 28 0.036 0.003 ± 0.003 128 3 

Al - -0.04 30 0.037 0.064 ± 0.020 128 3 

A2 - -0.08 31 0.039 0.110 ±0.016 128 3 

A3 - -0.14 33 0.041 0.152 ±0.040 128 3 

A4 -0.17 37 0.047 0.172 ±0.040 128 3 

~B1 -50.00 -0.01 -0.19 33 0.042 0.149 ±0.022 \W 

B2 -25.00 -0.02 -0.19 34 0.042 0.151 ±0.014 128 3 

B3 -10.00 -0.04 -0.20 32 0.040 0.176 ±0.017 128 3 

B4 -5.00 -0.09 -0.22 29 0.036 0.147 ±0.012 128 3 

B5 -2.50 -0.18 -0.23 28 0.035 0.1 19 ±0.021 128 3 

B6 -1.00 -0.49 -0.25 26 0.032 0.069 ±0.010 128 3 

B7 -0.50 -1.02 -0.25 25 0.031 0.033 ±0.001 128 3 

B8 -0.25 -2.25 -0.28 23 0.028 0.014 ±0.002 128 3 

B9 -0.10 -7.67 -0.38 17 0.021 0.015 ±0.005 128 3 

CI L99 O20 -0.20 32 0.040 0.051 ±0.013 128 3 " 

C2 1.75 0.24 -0.21 30 0.038 0.039 ± 0.024 128 3 

C3 1.50 0.30 -0.22 28 0.036 0.053 ± 0.003 128 3 

C4 1.25 0.38 -0.24 26 0.033 0.076 ±0.014 128 3 

C5 1.00 0.49 -0.25 26 0.032 0.066 ± 0.006 128 3 

C6 0.75 0.67 -0.25 25 0.032 0.047 ± 0.007 128 3 

C7 0.50 1.02 -0.25 25 0.031 0.014 ± 0.006 128 3 

C8 0.25 2.19 -0.27 23 0.029 -0.015 ± 0.005 128 3 

C9 0.10 7.46 -0.37 17 0.021 -0.036 ±0.001 128 3 
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TABLE 2 

Summary of the Runs with Keplerian Shear, i.e. here q - 



1.5 in all Runs. 



Run 


Pr 


Ra 


Co 


Sh 


Re 


Ma 




grid 


Dl 


j 


10 


0.30 


-0.22 


28 


0.036 


0.053 ± 0.003 


128 3 


D2 




2 • 10 6 


0.29 


-0.21 


59 


0.037 


0.053 ± 0.028 


256 3 


D3 




4- 10 6 


0.27 


-0.20 


125 


0.039 


0.033 ± 0.004 


256 3 


D4 


To 


10 7 


0.27 


-0.20 


315 


0.039 


0.018 ±0.006 


512 3 


D5 


in 


2- 10 7 


0.26 


-0.20 


648 


0.041 


0.031 ±0.010 


1024 3 


El 




10 b 


0.01 


-0.01 


27 


0.034 


0.022 ± 0.021 


128 3 


E2 




10 6 


0.03 


-0.02 


26 


0.033 


0.032 ± 0.022 


128 3 


E3 




10 6 


0.06 


-0.05 


28 


0.035 


0.056 ± 0.013 


128 3 


E4 




10 6 


0.12 


-0.09 


27 


0.034 


0.084 ± 0.007 


128 3 


E5 




10 6 


0.29 


-0.21 


29 


0.037 


0.045 ±0.013 


128 3 


E6 




10 6 


0.63 


-0.47 


27 


0.034 


0.074 ±0.011 


128 3 


E7 




It) 6 


1.32 


-0.99 


25 


0.032 


0.041 ±0.011 


128 3 


E8 




10 6 


4.33 


-3.24 


20 


0.025 


-0.040 ± 0.003 


128 3 



gradient matrix of the large-scale velocity 
R ij = -MijkiUkj + ..., 



(17) 



where the dots indicate higher derivatives of U that can occur 
in the expansion. The expression dTTb states that the stress 
is diffusive in character. In the general case the fourth rank 
tensor A/y*/ can have a complicated structure, see e.g. Riidiger 
(1989). However, if we assume that the shear is weak, the 
simplest description of the horizontal stress generated by our 
linear shear flow is given by 



R 



-v t (U XJ + U y , x ) = -2v t S. v . v = -v t S, 



(18) 



where v, = v t (z) is the z-dependent turbulent viscosity, and 
the component S vv of the rate-of-strain tensor is not to be con- 
fused with the shear rate S . To our knowledge, SB96 present 
the only published results of turbulent viscosity in the ab- 
sence of rotation as determined from convection simulations 
with large-scale shear, and they only provide a volume av- 
eraged quantity for one case. If also rotation is present (cf. 
Cabot 1996; Lesur & Ogilvie 2010) the Reynolds stress can 
no longer be related to turbulent viscosity alone (see below). 
In the present paper we study the dependence between stress 
and shear systematically and estimate the turbulent viscosity 
coefficient v t . 

It turns out that in many applications, Equation (TTTb is in- 
sufficient to describe the stress. For example, according to 
Equation ([Tol l, the Reynolds stress component Rg^ derived 
from observations of sunspot proper motions with the ob- 
served surface differential rotation would yield v t < 0, which 
is clearly unphysical (see the discussions in Tuominen & 
Riidiger 1989; and Pulkkinen et al. 1993). This motivates the 
inclusion of a non-diffusive contribution proportional to the 
rotation of the system (e.g., Wasiutynski 1946), such that 



(19) 



where Ay* are the components of the A-effect. This effect 
is expected to occur in anisotropic turbulence under the influ- 
ence of rotation (e.g. Riidiger 1989). In convection the density 
stratification provides the anisotropy. This is confirmed by 
numerous simulations of rigidly rotating stratified convection 
(e.g. Pulkkinen et al. 1993; Chan 2001; Kapyla et al. 2004; 
Riidiger etal. 2005). Although additional shear flows are gen- 
erated in these systems when the gravity and rotation vectors 



are not parallel or antiparallel, no serious attempt has been 
made to quantify the turbulent viscosity in convection. 

When shear and rotation are both present, it is no longer 
possible to distinguish between the diffusive and non-diffusive 
contributions without resorting to additional theoretical argu- 
ments. Here we use a simple analytical approach based on 
the so-called minimal r-approximation (see e.g. Blackman & 
Field 2002, 2003) to estimate the contributions from v t and 
the A-effect. 

The idea behind the minimal r-approximation is to use re- 
laxation terms of the form -T~ x Rij in the evolution equations 
for the components of the Reynolds stress Rjj as a simple de- 
scription of the nonlinearities. Using the decomposition < fT3T > 
and the Navier-Stokes equations one can derive equations for 
the Reynolds stress. For the purposes of the present study it 
suffices to consider a situation with one spatial dimension, z. 
In this case the evolution equations can be written as 

d,R u = -V ii R :! - R iz d z Uj - Rj-dJJi -S6 yi R xj 

-S 5 yj R xi - 2eu k Q.iR kj - 2e jlk Q.,R ki + N u , (20) 

where TYy represents nonlinear terms. As described above, us- 
ing the minimal r-approximation as closure model, one sub- 
stitutes 

N u = -T-'Rij, (21) 

where r is a relaxation time, which is here treated as a free pa- 
rameter that we expect to be close to r tum . In the present paper, 
however, we will not solve the closure model self-consistently 
but rather compare the leading terms with the numerical re- 
sults. A more thorough investigation using the full closure 
model is postponed to a later study. 

4. RESULTS 

In order to study the Reynolds stress generated by shear 
and rotation, we perform five sets of simulations summarized 
in Tables [T]and |2] Our base run is one with Re w 30, Pr = 1, 
and Co = Sh = (run AO). In set A we add only shear and in 
sets B and C we keep the shear constant and vary rotation so 
that q is negative and positive, respectively. In the remaining 
sets of calculations we study the Keplerian shear case (q = 
1.5): in set D we vary Re and thus Ra with fixed shear and 
rotation whereas in set E we vary both Co and Sh with fixed 
Re. 

4.1. Only shear, Sh + 0, Co = 

When only shear is present, the turbulent viscosity can be 
computed from Equation ( fT8l ) as 

v t = -R^/S, (22) 

where the superscript S indicates that the stress is due to the 
large-scale shear. By virtue of density stratification, the com- 
ponents of the stress tensor, and hence v t , are functions of 
depth, i.e. 2?y = Rij(z) and v t = v t (z). We normalize our re- 
sults with the estimate 



VtO 



(23) 
(24) 
(25) 

Note that, if we allow St + 1, we have v t o = ^Stu lms k^ 1 and 
the ratio Vt/v t o gives a measure of the Strouhal number. A 



Assuming that the Strouhal number 

St = TUnnskf, 

is equal to unity, i.e. r = rtum, we obtain 



Vto - ^u ms k { . 
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Fig. 2. — Vertical pro file o f the turbulent viscosity v, (thick solid line), 
computed from Equation ^22t. and the analytical estimate (thick dashed line) 
according to Equation j28t for run Al with Sh = -0.04, Co = and Re = 
30. The shaded areas show the error estimates. The dotted vertical lines at 
z = and z = d denote the bottom and top of the convectively unstable layer, 
respectively. 
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Fig. 3. — Turbulent viscosity v t as a function of height for different shear pa- 
rameters (runs A1-A4). The inset shows R„ averaged over the convectively 
unstable layer as a function of Sh; dotted line is proportional to Sh. 



typical example of v t is shown in Figure [2] We find that the 
stress owing to the large-scale shear is always positive, yield- 
ing v t > 0. The stress R]^ increases roughly proportionally to 
the shear parameter S in the range 0.04 < |Sh| < 0.17 so that 
the ratio v t /v t o remains approximately constant; see Figure[3] 
Increasing the shear much beyond our relatively low maxi- 
mum value leads to large-scale vorticity generation. These 
flows often grow so strong that the Mach number reaches 
unity leading to numerical difficulties. We associate this phe- 
nomenon with the 'vorticity dynamo' observed in forced tur- 
bulence simulations (Yousef et al. 2008; Kapyla et al. 2009a) 
and several theoretical studies (e.g. Elperin et al. 2003, 2007). 
Although the large-scale vorticity generation is weak in our 
runs A1-A4 we see that in the absence of rotation even mod- 
est shear increases the rms-velocity measurably (see Table[T|i. 
However, when rotation is added, this instability vanishes 
(see, e.g., Snellman et al. 2009). 

Our estimate for v t o is based on the volume averaged rms- 
velocity and a somewhat arbitrarily defined length scale 1 /kf . 
These choices are partly responsible for the large values of 
v t /vto- In order to obtain a more accurate estimate we derive 



an evolution equation for R xv using Equation ( f20b 

d,R xy = -U z d z R xy - RyzdJJx - R.xzd-Uy - zt|S + Nxy. (26) 

In general the first term on the rhs is non-zero in the com- 
pressible case but, as our Mach numbers are small, U z is neg- 
ligible. Furthermore, the stress components R xz and R yz van- 
ish under the assumption that the imposed shear is the only 
large-scale velocity component that depends on the horizon- 
tal coordinates. Thus the only terms remaining are 

d,R xy = -4S + N xy . (27) 

We now ignore the time derivative and apply the minimal r- 
approximation, i.e. Equation ( 1211 ). to obtain 

R xy = -tuIS = -v t S. (28) 

Note that we use t as a fitting parameter when comparing the 
different sides of the equation. We find that the simple an- 
alytical estimate can be fitted with the stress from the simu- 
lations when the Strouhal number is in the range 3 ... 4 for 
all runs in set A, e.g., see the comparison shown in Figure [2] 
for run Al. Here, for simplicity, we have assumed that t has 
no dependence on z which can contribute to the fact that the 
curves have somewhat different depth dependencies. Strouhal 
numbers in the range 1-3 are in line with previous numerical 
findings with forced turbulence (e.g. Brandenburg et al. 2004) 
and convection (e.g. Kapyla et al. 2009b). 

4.2. Shear and rotation, Sh, Co + 
Figures|4]and|5]show the total stress R„ from sets B and C 
where the imposed shear with S = -0.05 -\/g/d is kept con- 
stant and rotation is varied in a way that q is either negative 
(set B) or positive (set C), respectively. For shear parameters 
q > 2, the flow is Rayleigh unstable; thus, we investigate the 
parameter regime from slightly below 2 down to larger nega- 
tive values. Note that although the imposed shear is constant, 
the value of Sh varies somewhat because it is based on the 
turbulent velocity u rms which itself is a function of shear and 
rotation. 

We find that for negative q (Figure HJ, the stress decreases 
monotonically as rotation is increased. For slow rotation, i.e. 
Co > -0.04, the differences between the runs are not sta- 
tistically significant, cf. Figure [2] and Table [T] This is also 
clear from Figure [6] where we plot the volume average of the 
stress over the convectively unstable layer as a function of ro- 
tation. For more rapid rotation, we interpret the decrease of 
the stress as the non-diffusive contribution due to the A-effect 
which is likely a good approximation for slow rotation (see 
below). The situation is less clear in the case of positive q, 
see Figure|5]and the inset of Figure|6] A contributing factor is 
that we cannot use arbitrarily small Q. in the positive q regime 
because cases with q > 2 are Rayleigh unstable. However, 
in the runs that can be done the stress is somewhat decreased 
in comparison to runs with corresponding |Co| and negative 
q. For the most rapidly rotating cases in set C the sign of 
the stress changes, which is not observed in set B. The sign 
change suggests that the A-effect dominates over the turbu- 
lent viscosity in the rapid rotation regime. Similar asymmetry 
between the regimes of positive and negative q have also been 
reported e.g. by Snellman et al. (2009) and Korpi et al. (2010). 

It is interesting to see whether simple analytical models can 
reproduce the simulation results. For example, in the presence 
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Fig. 4. — Total stress R xy for different rotation rates using a subset of runs 
in set B with negative q, Sh as -0.22, and Re = 30. 
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Fig. 5. — Same as Figure|4] but for a subset of runs from set C with positive 
values of q (see the legend). 
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Fig. 6. — Total stress R„, averaged over the convectively unstable layer for 
the runs in set B. The dashed and dot-das hed curves show results from the 
minimal r approximation model (Equation [29}. The open symbol on the left 
denotes the stress for run A4 with no rotation. The inset shows the same 
representation for the runs in set C. 



(R„) = a x St 



-(Co + Sh) (R ( ®) + Co(Rfh 



(29) 



of weak shear and rotation, the minimal r-approximation ap- 
plied to a homogeneous system with no convection gives (see 



1 +4CoSt 2 (Co + Sh) 

where angular brackets denote volume averaging and the su- 
perscript refers to a non-rotating and non-shearing reference 
state. A similar result was obtained earlier for arbitrary shear 
and rotation by Narayan et al. (1994) with a conceptually dif- 
ferent model where individual eddies were treated as parti- 
cles that scatter off each other. We note that neither of these 
models is directly applicable to the present system, although 
treatment of convection can be introduced into the models 
(e.g. Kumar et al. 1995). However, our purpose here is not 
to perform a detailed comparison of the closure models with 
simulations but rather to attain a broad understanding of the 
system. In Figure [6] we compare the numerical results with 
the analytical estimate given in Equation ( 1291 . keeping the 
Strouhal number as a free parameter. Furthermore, we have 
introduced a scaling parameter a\ (=1.7) in order to improve 
the fit. 

We find that parameters St = 1.0 and a\ = 1.7 produce 
a good fit to the numerical results for the runs in set B. We 
have here normalized our results from Equation ( 1291 by the 
value of M rms from the non-rotating run AO, so a scaling factor 
ct\ greater than unity can in principle be understood to reflect 
the decrease of u^ as rotation is increased. However, this 
scaling is not essential since even the unsealed curve shows 
qualitatively the same behavior. The fit for the runs in set C 
is not as successful although the simple model coincides with 
the simulation data for intermediate values of Co. However, 
the simulation results fall below the model for q > 1 .25 and 
the negative stresses for rapid rotation are not captured by the 
model. The latter is in line with the discovery of Snellman 
et al. (2009) that the validity of the minimal r-approximation 
is limited to the slow rotation regime. The lack of proper pa- 
rameterization of thermal convection in our simple model is 
another obvious reason for the differences. 

4.3. A-effect due to shear-induced anisotropy 

In the absence of shear, but including rotation parallel or 
antiparallel to gravity, turbulence is statistically axisymmetric 
and there is no asymmetry between the turbulence intensities 
in the two horizontal directions, i.e. R xx = R yy . This implies 
that there is no horizontal A-effect which, to the lowest order, 
is proportional to (Riidiger 1989) 

Rf = A H Q = 2Q () r(/?, v - R xx ). (30) 

Note that the same result is borne out of Equation d20b if we 
assume that all large-scale flows vanish and allow deviations 
from axisymmetry, i.e. R xx + R yy . However, when shear is 
included, the turbulence becomes anisotropic in the horizon- 
tal plane, enabling the generation of a non-diffusive contri- 
bution to the stress R xy due to rotation, according to Equa- 
tion d30b . Such contributions have earlier been studied ana- 
lytically (Leprovost & Kim 2007, 2008a,b) and numerically 
(Snellman et al. 2009) for isotropically forced homogeneous 
turbulence. 

When both shear and rotation are present it is not possible 
to separate the diffusive from the non-diffusive contribution 
directly. Furthermore, using the diffusive stress from a purely 
shearing run to extract the non-diffusive one from a run with 
both the effects is also problematic due to the relatively large 
errors in the data (cf. Figure [6]) which can lead to spurious 
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results. The large errors in the purely shearing runs can pos- 
sibly be explained by a subcritical vorticity dynamo (cf. Sec- 
tion l4.ll ). However, if rotation is slow, we can use the simple 
analytical results of Equations d28l l and d30T > to express the 
total stress as 



R xy = 2D.T(Ryy - R xx ) - tR xx S. 



(31) 



On the other hand, we can express the stress in terms of the 
A-effect and turbulent viscosity by 



Rxy — AhO - v t S , 



(32) 



where we now have 



A H = 2r(Ryy - R xx ), v, = tR x 



(33) 

and where we again treat t as a fitting parameter. Figure [7] 
shows an example from run B5. We find that, when r corre- 
sponding to St 2 is used, the total stress is in broad agree- 
ment with Equation (13 It . In addition to the weaker negative 
diffusive contribution, corresponding to the turbulent viscos- 
ity, we find a non-diffusive part of the opposite sign. Reason- 
ably good fits can be obtained for runs with |Co| < 1 with a 
r corresponding to St « 2. For more rapidly rotating cases 
the representation Equation (l3TT l breaks down. Furthermore, 
in the rapid rotation regime the relevant time scale is the rota- 
tion period rather than the turnover time. For the purposes of 
visualization, and without altering the qualitative character of 
the results, we consider here volume averages over the con- 
vectively unstable layer. Results for runs B1-B7 are shown 
in Figure [8] We use a fixed r = 8 -\/d/g which corresponds 
to St «; 1.5 .. .2.1 in these runs. We find that the total stress 
is fairly well reproduced for runs where |Co| is below unity. 
Furthermore, the diffusive contribution stays almost constant 
as a function of Co, whereas the non-diffusive part is close to 
zero for |Co| < 0.1. The turbulent viscosity, as obtained from 
Equation d33t . shows a weak declining trend as a function of 
rotation, see the lower panel of Figure [8] The coefficient Ah 
has values in the range (0.5-l)v t o for slow rotation. The error 
estimates increase towards slow rotation which is consistent 
with the fact that the non-diffusive stress is small at low Cori- 
olis numbers. These results suggest that the A-effect is non- 
zero when the anisotropy of the turbulence is induced by the 
shear flow with a roughly rotation independent coefficients 
Ah- However, as our method breaks down when |Co| > 1, 
we cannot study quenching behavior of v t and the A-effect for 
rapid rotation. 

4.4. Dependence on Reynolds number 

In set D we vary the Reynolds number, keeping shear and 
rotation fixed. Here we choose q — 1.5, corresponding to 
Keplerian shear. Again, the values of Co and Sh are not ex- 
actly constant due to the varying u lms ; see Table [2] Figure [9] 
shows the results for the horizontally averaged stress from the 
runs in set D. We find that for relatively weak shear and rota- 
tion, the stress is positive in the convectively unstable layer — 
in apparent contradiction with some earlier results (Stone & 
Balbus 1996; Cabot 1996) but in accordance with the recent 
results of Lesur & Ogilvie (2010). We discuss this issue in the 
next section in detail. 

We also find that the vertical distribution and magnitude of 
the stress are not much affected when the Reynolds number 
is increased from 28 to 648, see the inset of Figure [9] where 
Rxy averaged over the convectively unstable layer and time is 
shown. We note that in addition to the Reynolds number, the 
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Rayleigh number in this set changes by a factor of 20. How- 
ever, we have kept the heat conduction, K, and thus the input 
energy flux, constant so that M rms varies by only 10% within 
set D. Had we kept Pr = 1, the energy input at the lower 
boundary would have also decreased by a factor of 20. This 
would have resulted in a much lower w lms and thus propor- 
tionally greater values of Co and Sh. This would have likely 
produced a very different trend as a function of Ra because Co 
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Fig. 9. — Horizontally averaged stress component R„, from set D with 
varying Reynolds number and Co » 0.3 and Sh = -0.2. The inset shows 
the stress averaged over the convectively unstable layer as a function of the 
Reynolds number. 
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Fig. 10. — Vertical profiles of S xy for several values of Co using data for 
set E. The inset shows R xy averaged over the convectively unstable layer as a 
function of Co. 



and Sh are considered as the relevant dimensionless parame- 
ters for the Reynolds stress (see, e.g., Snellman et al. 2009). 

In a recent study, Lesur & Ogilvie (2010) presented results 
from Boussinesq convection in an otherwise similar shearing 
box setup as ours. They report that the stress changes sign 
from negative to positive in the Keplerian case for strong shear 
when the Rayleigh number is increased sufficiently. This ap- 
pears to be in contradiction with our results regarding the de- 
pendence on Rayleigh number. However, in their setup the 
Richardson number (Ri), defined as the negative of the ratio 
of the squared Brunt-Vaisala frequency and the squared shear 
rate, is less than unity. This indicates that their models are in 
the shear-dominated regime whereas in our case only one sim- 
ulation (run E8) falls in this regime. The low-Ri case at high 
Rayleigh numbers certainly merits further study. 

4.5. Relation to accretion disk theory 

In their paper, Stone & Balbus (1996) performed a numer- 
ical simulation of convection in a local accretion disk model 
and found that the total stress is small and on average directed 
inward. This numerical result based on one simulation and in- 
sight drawn from analytical arguments led them to conclude 
that convection cannot account for the outward angular mo- 
mentum transport required in astrophysical accretion disks. 



However, numerical simulations of isotropically forced ho- 
mogeneous turbulence under the influence of shear and rota- 
tion indicate that the total stress can change sign as a function 
of Co when q is fixed (Snellman et al. 2009). In their study, 
Snellman et al. (2009) found that, for slow rotation, the stress 
is positive and changes sign near Co = 1 . For rapid rotation 
(Co « 10) the stress appears to drop close to zero. In the con- 
text of mean-field hydrodynamics, this can be understood as 
quenching of the A-effect and turbulent viscosity due to shear 
and rotation, i.e. Ah = Ah(£2, S) and v t = v t (£2, S). 

In an effort to study whether these results carry over to con- 
vection we perform a set of runs where we keep q - 1.5 fixed 
and vary the values of Co and Sh. Our results for the hor- 
izontally averaged stress are shown in Figure [10] For slow 
rotation and weak shear R„, is positive with a maximum value 
of roughly 10 per cent of the mean square velocity. As we in- 
crease the rotation, the stress decreases and changes sign for 
Co * 2. For Co = 4.33, R xy is negative everywhere in the 
convectively unstable layer. This is when the flow pattern has 
become markedly anisotropic, as can be seen from visualiza- 
tions of U z on (or near) the periphery of the computational 
domain (Figure ITTb. Indeed, the flow pattern becomes rather 
narrow in the cross-stream direction, while being roughly un- 
changed in the streamwise direction. In the rapid rotation 
regime the non-axisymmetric structures tend to disappear and 
promote a two-dimensional flow structure that leads to inward 
transport (Cabot 1996), which is also visible from our results. 

Let us now ask under what conditions can we expect out- 
ward angular momentum transport in accretion disks, i.e., 
when can one expect the Coriolis number to be below unity? 
In disks we have CloH = c s , where H is the scale height and 
c s is the sound speed. Inserting this into our definition of Co 
we obtain from Equation ( fT2l 

2Q 2O H 1 2 i 
Co = = — = Ma t , (34) 

Wrms^f Mrms kfH kfH 

where we have defined the turbulent Mach number Ma t = 
Wrms/cs- In our simulations, kfH is of order unity, which would 
suggest that we can expect Co < 1 only for supersonic turbu- 
lence. 

A more refined estimate can be obtained by invoking the 
Shakura-Sunyaev a-parameter (Shakura & Sunyaev 1973), 
which is introduced as a parameterization of the turbulent vis- 
cosity v t via 

v t = ac s H. (35) 

We can use this parameterization together with Equation d23l 
to eliminate M rms in favor of a under the assumption that v t ~ 
v t rj. This leads to 

2c s 2 
C ° " 3^H ~ 3a(W (36) 

In our simulations we have H » 0.62a! at the base of the con- 
vectively unstable layer. Together with our estimate kf -2n/d 
this yields kfH « 4, and therefore Co w (24a)~ 1 . We can 
therefore expect Co < 1 for accretion disks where the turbu- 
lence is sufficiently vigorous so that a > 0.04. 

We now use our simulations to estimate a. In disks, the rate 
of strain is proportional to qClo, so the total turbulent stress is 
given by 

T v = vtfOo, (37) 

where, in the absence of any other stresses such as from mag- 
netic fields, the total stress per unit mass is given by 

Txy = Rxy = U x U y . (38) 



Fig. 11. — Effect of increasing rotation and shear rates on the visual appearance of the vertical velocity U L (runs E2, E4, E6 and E8). All runs are for Keplerian 
shear, i.e. q = 3/2 and hence Sh/Co = -3/4, the resolution is 128 3 mesh points and the Reynolds number varies between 20 and 30. 
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Fig. 12. — Viscosity parameter a according to Equation i39\ for the runs 
in set E. The inset show the same in log-log scales. The dashed line shows a 
scaling inversely proportional to the Coriolis number for reference. 



Equation d36l ). We find that for Co « 1 we have a ~ 4 • 10 
which is two orders of magnitude smaller than the estimate 
derived above. 

Another problem facing the suggestion that convection 
might drive the angular momentum transport in accretion 
disks is that without an internal heat source in the disk, con- 
vection is not self-sustained (cf. SB96). However, many disks 
are likely to be susceptible to the magnetorotational instability 
which can extract energy from the shear flow and ultimately 
deposit it as thermal energy in the disk. If the material in 
the disk is sufficiently thick optically, the pile-up of energy 
from the magnetorotational instability could render the verti- 
cal stratification of the disk convectively unstable. 

As alluded to in the introduction, our aim is not to claim 
that convection is solely responsible for the outward angular 
momentum transport in accretion disks but to show that, given 
the right conditions, convection can contribute to outward an- 
gular momentum transport. 



Combining Equations d35l l and (l37l i gives 

T„ 



(39) 



We compute T xy and hence a as volume averages over the 
convectively unstable region. For the normalization factor, 
we take conservative estimates of c s and H from the bottom 
of the convectively unstable layer. We find that for slow rota- 
tion, i.e. Co < 0.2, a is roughly constant with a value of the 
order of 0.01 (see Figure fl2l). For more rapid rotation a de- 
creases and eventually changes sign. The points in the range 
Co w 0.06 ... 1 are consistent with a scaling inversely pro- 
portional to the Coriolis number which is also suggested by 



5. CONCLUSIONS 

The present results have shown that hydrodynamic convec- 
tion is able to transport angular momentum both inward and 
outward, depending essentially on the value of the Coriolis 
number, in accordance with earlier results from homogeneous 
isotropically forced turbulence (Snellman et al. 2009). This 
underlines the importance of considering comprehensive pa- 
rameter surveys and not to rely on demonstrative results from 
individual case studies. For given value of the Coriolis num- 
ber, the stress is found to be relatively independent of the 
value of the Rayleigh number (Section l4~4b . By varying shear 
and rotation rates separately, we have been able to quantify 
the relative importance of diffusive and non-diffusive contri- 
butions to the Reynolds stress tensor. In agreement with ear- 
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lier work, it turns out that the turbulent kinematic viscosity is 
of the order of the mixing length estimate and has roughly the 
same value as the turbulent magnetic diffusivity found earlier 
for similar runs (Kapyla' et al. 2009b). In other words, the 
turbulent magnetic Prandtl number is around unity, again in 
accordance with results from simpler systems (e.g. Yousef et 
al. 2003). 

The other important turbulent transport mechanism in rotat- 
ing turbulent bodies is the A effect. Although the importance 
of this effect is well recognized in solar and stellar physics 
(e.g. Riidiger & Hollerbach 2004), it is not normally consid- 
ered in the context of accretion disks. In the present paper we 
have been able to quantify its importance for a range of Cori- 
olis numbers by means of a simple analytical model making 
use of the minimal r-approximation. For slow rotation the 
coefficient A H is of the order of v t0 and independent of the 



Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214 
Blackman, E. G. & Field, G. B. 2002, PhRvL, 89, 265007 
Blackman, E. G. & Field, G. B., 2003 Phys. Fluids, 15, L73 
Brandenburg, A., Kapyla, P.J. & Mohammed, A. 2004, Phys. of Fluids, 16, 
1020 

Cabot, W. & Pollack, J.B. 1992, GAFD, 64, 97 
Cabot, W. 1996, ApJ, 465, 874 
Chan, K. L. 2001, ApJ, 548, 1102 

Elperin, T., Kleeorin, N. & Rogachevskii, I. 2003, PRE, 68, 016311 
Elperin, T., Golubev, N., Kleeorin, N. & Rogachevskii, I. 2007, PRE, 76, 
066310 

Korpi, M.J., Kapyla, P.J. & Vaisalii, M.S., 2010, AN, 331, 34 
Kapyla, P.J. & Brandenburg, A., 2008, A&A, 488, 9 
Kapyla, P.J., Korpi, M.J. & Tuominen, I., 2004, A&A, 422, 793 
Kapyla, P.J., Korpi, M.J. & Brandenburg, A., 2008, A&A, 491, 353 
Kapyla, P.J., Korpi, M.J. & Brandenburg, A., 2009b, A&A, 500, 633 
Kapyla, P.J., Mitra, D. & Brandenburg, A., 2009a, PRE, 79, 016302 
Krause, F. & Riidiger, G, 1974, AN, 295, 93 

Krause, F. & Radler, K.-H., 1980, Mean-field Magnetohydrodynamics and 

dynamo theory (Pergamon Press, Oxford), 27 1 
Kumar, P., Narayan, R. & Loeb, A., 1995, ApJ, 453, 480 
Narayan, R., Loeb, A. & Kumar, P., 1994, ApJ, 431, 359 
Leprovost, N. & Kim, E.J. 2007, A&A, 463, 9L 
Leprovost, N. & Kim, E.J. 2008a, PRE, 78, 016301 



Coriolis number. However, once the Coriolis number exceeds 
a value around unity, our method of separating the turbulent 
viscosity and the A-effect breaks down, which reinforces the 
need for a truly independent determination not only of diffu- 
sive and nondiffusive contributions to the Reynolds stress, but 
also of all the components of the full stress tensor. 



The computations were performed on the facilities hosted 
by CSC - IT Center for Science Ltd. in Espoo, Finland, 
who are administered by the Finnish Ministry of Education. 
Financial support from the Academy of Finland grants No. 
136XYZ, 140970 (PJK) and 112020 (MJK) are acknowl- 
edged. The authors acknowledge the hospitality of NORDITA 
during their visits. 



Leprovost, N. & Kim, E.J. 2008b, PRE, 78, 036319 

Lesur, G. & Ogilvie, G.I. 2010, MNRAS, 404, L64 

Riidiger, G, Tschape, R. & Kitchatinov, L.L., 2002, MNRAS, 332, 435 

Pulkkinen, P., Tuominen, I., Brandenburg, A., Nordlund, A & Stein, R. F., 

1993, A&A, 267, 265 
Riidiger, G. 1989, Differential Rotation and Stellar Convection: Sun and 

solar-type stars (Berlin, Akademie-Verlag) 
Riidiger, G, & Hollerbach, R., 2004 The Magnetic Universe (Wiley-VCH, 

Weinheim), 47 
Riidiger, G., Egorov, P. & Ziegler, U. 2005, AN, 326, 315 
Ryu, D. & Goodman, J., 1992, ApJ, 388, 438 
Shakura, N.I. & Sunyaev, R.A., 1973, A&A, 24, 337 

Snellman, J.E., Kapyla P.J., Korpi, M,J. & Liljestrom, A.J., 2009, A&A, 505, 
955 

Stone, J.M. & Balbus, S.A., 1996, ApJ, 464, 364 
Tuominen, I. & Riidiger, G. 1989, A&A, 217, 217 
Wasiutynski, J. 1946, ApNr, 4, 1 
Wisdom, J. & Tremaine, S. 1988, AJ, 95, 925 

Yousef, T. A., Brandenburg, A., & Riidiger, G. 2003, A&A, 411, 321 
Yousef, T. A., Heinemann, T., Schekochihin, A. A., et al. 2008, PhRvL, 100, 
184501 



Id: paper.tex.v 1.143 2010-07-02 08:27:08 pkapyla Exp 



